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ABSTRACT 



CN ! This paper studies multistep methods for the integration of reversible dy- 

namical systems, with particular emphasis on the planar Kepler problem. It has 
K* I previously been shown by Cano & Sanz-Serna that reversible linear multisteps 

f^^ I for first-order differential equations are generally unstable. Here, we report on a 

^ I subset of these methods - the zero-growth methods - that evade these instabili- 

O ! ties- We provide an algorithm for identifying these rare methods. We find and 

q} I study all zero-growth, reversible multisteps with six or fewer steps. This select 

"Th I group includes two well-known second-order multisteps (the trapezoidal and ex- 

Q-i| plicit midpoint methods), as well as three new fourth-order multisteps - one of 

O ' which is explicit. Variable timesteps can be readily implemented without spoil- 

^ I ing the reversibility. Tests on Keplerian orbits show that these new reversible 

multisteps work well on orbits with low or moderate eccentricity, although at 
least 100 steps/radian are required for stability. 



1. Introduction 

The success of symplectic integration algorithms (siAs) as tools for the numerical solu- 
tion of Hamiltonian systems illustrates the importance of numerical algorithms that inherit 
the fundamental physical constraints of the systems to which they are applied ( "geometric 
integrators" ) . 

Time- reversal symmetry plays a central role in physics (e.g., Davies 1974). To define 
reversibility formally and without reference to coordinates, let x denote the state of the 
system governed by the differential equations 

l^/W. (1) 
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Let T be an involution, i.e. an operator such that T^x = x. The system (|l]) is T-reversible 
if T reverses the direction of time, that is, if 

^ = -f(Tx). (2) 

In standard phase space x = (q, p) and T = diag (1,-1). 

The trajectory of the system governed by the differential equations (|l|) is defined by a 
map Gt such that a system located at x at time is found at GfX at time t; by definition 
Gt = Gz\. Then, T-reversibility implies that 

GtTGt = T or TGt = G^^T = G^tT. (3) 

Reversible systems may or may not be Hamiltonian. However, the general features of 
motion in reversible and Hamiltonian systems show many similarities, including the existence 
of families of KAM tori (e.g., Moser 1973; Arnold 1984; Roberts & Quispel 1992). These 
considerations motivate us to examine reversible integration algorithms (rias). A one-step 
numerical integration algorithm with timestep h is defined by an operator Gh that is intended 
to be a close approximation to Gh, and the algorithm is an RIA if when it is applied to a 
reversible system 

GhTGh = T. (4) 

Note that in contrast to Gt, Gh is not necessarily equal to GZh- 



1.1. Symplectic integration algorithms 

Let us first briefly review SlAs; for more detail see Channell & Scovel (1990), Yoshida 
(1993), Marsden et al. (1996), and Sanz-Serna & Calvo (1994). The standard example of 
a Hamiltonian system is motion in a conservative potential (this is also reversible). The 
Hamiltonian is 

H{q,p) = lp' + U{q) (5) 

and the equations of motion are 

^^ ^P T7TT (a\ 

or 
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The most popular SIA for such systems is the leapfrog or Verlet algorithm (e.g., Hockney 
& Eastwood 1981). This is a map from the phase-space coordinates (qn,Pn) at time t„ to 
coordinates (qn+i, Pn+i) at time tn+i = tn + h, defined by 

q' = qn. + ^hpn, Pn+1 = Pn " /iVf/(q'), qn+1 = q + |/iPn+l- (8) 

Leapfrog is an explicit second-order method (i.e. the one-step error is 0(/i^)), but higher- 
order explicit SlAs can be constructed by concatenating leapfrog steps of different sizes, 
including backwards steps. Leapfrog can be generalized to any separable Hamiltonian of the 
form H = J2Hj where each Hj is integrable. 

For general Hamiltonians, SlAs can be constructed by two methods. The first is based 
on expanding the scalar generating functions for the symplectic transformation (q„, p„) -^ 
(q„+i,p„+i) in powers of h (Kang 1986; Channell & Scovel 1990). The simplest example is 
the first-order SIA 

qn+1 = qn + h , p„+i = p„ - h , (9) 

which is implicit in general but explicit for Hamiltonians of the form (^). The second 
approach is to construct symplectic implicit Runge-Kutta algorithms. For the differential 
equation (|l|), the s-stage Runge-Kutta method is defined by (e.g., Ralston & Rabinowitz 
1978) 

s s 

ki = f{xn + hY^ aijkj), x„+i = x„ + /i ^ bjkj. (10) 

j=i i=i 

Runge-Kutta methods are symplectic if they satisfy a simple algebraic test (Sanz-Serna & 
Calvo 1994) 

biaij + bjaji-bibj = 0, l<i,j<s. (11) 

On setting i = j in (0), it is evident that all symplectic Runge-Kutta methods are necessarily 
implicit. The best known examples are the Gauss- Legendre Runge-Kutta methods (e.g., 
Sanz-Serna & Calvo 1994). The simplest can be written 

Xn+l = Xn -\- hj 2\-^n + Xn+l) , (12) 

which is the second-order implicit midpoint method. A symplectic fourth-order two-stage 
Runge-Kutta method is given by 

On = 022 = I, ai2 = i - IVS, a2i = 1 + ^Vs, 61 = 62 = |- (13) 

More generally, the Gauss- Legendre Runge-Kutta method is the unique s-stage method with 
order 2s, and this method is always symplectic. These attractive features are offset by the 
computational cost of solving the implicit equations ([T0|). 
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A limitation of SlAs is that they are difficuh to generahze to variable timesteps. When 
a standard variable timestep prescription is applied to an SIA, its performance is no better 
than that of non-symplectic integrators. The reason is simply that the mapping Gh{x){x) is 
usually not symplectic even when Gh{x) is. This is a serious problem, since most applications 
benefit from a variable timestep. One way to introduce variable timestep is by extending 
the phase space (Mikkola 1997). Suppose we wish to use a timestep g{c[,p). We define an 
extended phase space (Q, P) = ((go, q), {po, p)) by go = ^ and po = ~E, where t is time and 
E is energy. We then define a new Hamiltonian 

r(Q,P) = (7(q,p)[i/(q,p)+Po]. (14) 

The equations of motion for the Hamiltonian F in the extended phase space with fictitious 
time r are the same as the equations of motion for the Hamiltonian H in the original 
phase space, supplemented by the condition dt/dr = 5'(q, p). We may now integrate the 
Hamiltonian F using an SIA with fixed timestep At = 1, which corresponds to At ~ g. 
A limitation of this approach is that the Hamiltonian F is generally not separable, so that 
leapfrog and its generalizations cannot be applied. 



1.2. Reversible integration algorithms 

Several well-known one-step second-order algorithms for the differential equation (|l]) are 
RiAs: for example, the implicit midpoint method (eq. [l^), the trapezoidal method, 

Xn+l =Xn + \h [f{Xn) + f{Xn+l)] , (15) 

and the explicit midpoint method, 

Xn+2 = Xn + hf (x„+i) . (16) 

Some SlAs are also RiAs when they are applied to reversible Hamiltonians. Leapfrog and its 
generalizations are reversible. A reversible version of leapfrog with variable timestep g{q^, p) 
is given by (Huang & Leimkuhler 1997; Calvo et al. 1998) 

q' = qn + ;^Pn, P' = Pn. + l^Vf/(q'), 

2 h 

Pn+1 = , , ,. - Pn, Pn+1 = P' + 7^ Vf/(q'), 

^(q , P ) 2p„+i 

qn+1 = q' + 7, Pn+1, tn+1 = ^n + 7T (pn^ + Pnll) ■ (1^) 

^Pn+1 ^ ^ ^ 
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Quinn et al. (1997) describe a closely related algorithm based on a discrete set of timesteps 
that depend only on p or q and are separated by factors of two. 

A different approach to constructing RiAs is to modify non-reversible integration algo- 
rithms. For example, the following two operators are reversible even if Gh is not: 

G^sTG^/s^, {l + TGhT)-\l + Ghy, (18) 

if we write G/^ = 1 + hA + 0(/i^) then both of these operators are equal to 1 + ^h(A — TAT) 
to first order in h. Hut et al. (1997) describe numerical experiments with the second of these 
operators, calling it a "time-symmetrization meta-algorithm" since it can be applied to any 
one-step numerical integration algorithm. Reversible algorithms remain reversible with a 
variable timestep so long as the timestep is determined symmetrically by the location of the 
system at the start and end of the step, e.g. h = ^[g{xn) + g{xn+i)] (Hut et al. 1995). 

Special second-order differential equations of the form 

H T 

^-F(.) (19) 

are reversible, and so can profitably be integrated using RiAs. A useful source of high-order 
RiAs is linear multistep methods (e.g., Henrici 1962; Gear 1971; Lambert 1973), which have 
the form 

k 

J2 {ak-jXn^j - h%^jFj) = 0, (20) 

j=0 

where Xj = x(to + j^), Fj = F{xj). We may assume without loss of generality that a^ = \. 
The method is explicit if 6^ = and otherwise implicit. Linear multisteps include the classic 
Stormer (explicit) and Cowell (implicit) methods, which are characterized by at = 1, ak-i = 
—2, afc_2 = 1; and a/c_j = for j > 2. It is easy to show that the requirement for reversibility 
is that ttj = cttk-j and hj = cbk-j where c = ±1; thus the Stormer and Cowell methods 
are generally not reversible. Multistep RiAs are discussed by Lambert & Watson (1976), 
Quinlan & Tremaine (1990), and Fukushima (1998, 1999). In general, their performance 
over long time intervals on reversible dynamical systems is much better than Stormer- Cowell 
methods, although for occasional unfortunate choices of timestep their performance is ruined 
by timestep resonances (Quinlan 1999). 

Systems of first-order differential equations such as (|I|) can be integrated by linear 
multistep methods of the form 

k 
Y^ (ak-jXn-j - hPk-jfn-j) = 0, (21) 

3=0 
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where Xj = x(to + jh), fj = f{xj). We can assume without loss of generality that a^ = 1; 
explicit methods have (3^ = 0. These include the classic Adams-Bashforth (explicit) and 
Adams-Moulton (implicit) methods, which have a^-i = —1 and dk-j = for j > 2 (e.g., 
Henrici 1962; Gear 1971). Multistep methods for first-order differential equations are more 
general than multistep methods for special second-order equations, since any second-order 
equation can be written as a set of first-order equations. Moreover, implementing variable 
timesteps is easy in first-order equations - if we wish to use a timestep g{x) we introduce a 
fictitious time r by the relation dt = g{x)dT, and equation (|I|) can be rewritten as 

^ = 9{x)f{x), (22) 

which can be integrated using unit timestep in r. 

The aim of this paper is to investigate linear multistep RIAs for first-order differential 
equations. An important earlier investigation is that of Cano & Sanz-Serna (1998), who 
found that such methods typically possess grave numerical instabilities. We review general 
linear multistep methods in §^, multistep RiAs in § |2.1| , together with their instabilities in 



We show in §y that it is possible to construct some linear multistep RiAs that are not 
subject to the Cano & Sanz-Serna instabilities. A general discussion of stable multistep RiAs 



with up to 6 steps is given in §3.1 - §p73|. Finally, §H describes numerical examples based on 



integrating Kepler orbits and §| discusses our results. 

2. Linear multistep methods 

Following Henrici (1962) and Lambert (1973), we associate with ( pT]) the linear operator 

k 

L[xit), h]=J2 {cyk-jx[t + {k- j)h] - h(3k-jx[t +{k- j)h]}. (23) 

i=o 

We can expand x{t) in a Taylor series to obtain 

oo 

L[x{t),h]=T.C,h''x^''Kt), (24) 

g=0 



1 fc -| fe 



where 

k 

C, = Y.ai, C, = ^^a,/^--^— ^A/''"\ g = l,2,... (25) 

The order p of the multistep is the unique integer for which Co = ■ ■ ■ = Cp = 0, Cp+i ^ 0. 
The constant C = Cp+i/cr(l) is known as the error constant and is a measure of the local 
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truncation error. Now define the characteristic polynomials 

fc k 

p(0 = E«.e\ ^(0 = E/5.e- (26) 

If x(t) = exp(At), then 

L[x{t), h] = e^'[p{e^'') - Xhaie^'')], (27) 

which together with equation (|^) implies 

p(l + z)- log(l + z)a{l + z) = Cp+izP+^ + OizP+^), (28) 

where p is the order and z = exp(A/i) — 1. Order > requires 

fc 
p(l) = ^a, = 0, (29) 

3=0 

and order > 1 requires 

p'{l)=a{l). (30) 

A multistep method is zero-stable if and only if the roots of p{C,) lie inside the unit circle in the 
complex plane, or are on the unit circle and simple (proofs of this are given by Henrici 1962, 
Gear 1971 and Lambert 1973). Zero-stability ensures that the parasitic solutions generated 
by the additional roots of the difference equation (which is of order k, while the original 
differential equation (|I]) has order one) do not grow, at least in the limit of zero timestep. 
We will discuss other forms of stability shortly in §p.2|. 



We assume that the polynomials p{C,) and a{^) have no common roots other than 1, for 
the following reason. Suppose that ^g 7^ 1 is a common root, so that {C, — C,o) is a common 
factor. Then p(^) = (^ — ^o)p(O) ^(0 = (C ~ Co)5"(0 where p{^) and cr(^) are polynomials 
of degree k — 1. Then equation (pSJ) may be written 

p(l + z)- log(l + z)ail + z) = ^^^P+i + 0(^^+2). (31) 

1 — ^0 

Thus p{^), 5"(,^) define a (A; — l)-step method with the same order and the same error constant 
(C = Cp+i/(T(l) = Cp+i/[(l— ^0)5^(1)]) as the original fc-step method, and there is no obvious 
reason why the simpler method should not be used instead. 

There are 2k + 1 unknown coefficients in equation (^) (since a^ = 1), and therefore in 
principle these coefficients can be chosen so that the order is 2k (or 2A; — 1 if the method is 
explicit, with /J^ = 0). However, the maximum order of a zero-stable multistep method is 
A; + 1 if A; is odd and A; + 2 if A; is even (Henrici 1962). 



2.1. Reversible multistep methods 

Suppose that a trajectory {x„_fc, x^-k+i, ■■■, Xn}, {fn-k, fn-k+i, ■ ■ ■ , fn} satisfies the 
difference equation (PTl). The same segment of the time- reversed trajectory is given by 
{Txn,Txn-i, ■ ■ ■ ,Txn-k}, {On, dn-i, ■ ■ ■ , Qn-k} , whcre T is the time-reversal operator and 
Qk = f{Txk). We shall assume that Xk = (rfc,Vfc) and that fk = (v^, — V$(rfc)). Then 
T = diag(l, —1) and f{Tx) = —Tf{x). In an RIA the time-reversed trajectory should also 
satisfy the difference equation (|2T|) , that is 

k 

^{ak-jTxn-k+j - hf3k~j9n-k+j) = 0. (32) 

j=0 

Since gj = f{Txj) = —Tfj, we may operate with T to obtain 

k 

2_^{0ik-jXn-k+j + hPk-jfn-k+j) = 0, (33) 

j=0 



or, equivalent ly, 



k 
J2ia,Xn-j + h(3,fn-i) = 0- (34) 

i=o 



If (|3^) is to be satisfied whenever (|2T|) is satisfied, then we must have 

ak-j = caj, (3k- j = -cf3j, (35) 

where c is a constant. Applying this relation twice gives c^ = lsoc = ±l. Ifc = +1, 
we shall say that the multistep method has even parity; if c = — 1, the parity is odd. The 
characteristic polynomials (^) now satisfy 

p(0 = c^'pir'), ^(0 = -ceV(r'). (36) 

Now let ^j,j = 1, ... ,k he the roots of p(0- Equation (|36|) implies that if $,j is a root then so 
is ^~^. For zero-stability, the roots cannot lie outside the unit circle - therefore they must lie 
on the unit circle, and moreover they must be simple. Equation (|29|) guarantees that ^i = 1 
is a root for any method with order > 0. This root is simple if p'{l) 7^ 0, and equation ( pO] ) 
then requires that cr(l) 7^ for any method with order > 1. However, even-parity methods 
have cr(l) = Yl'j=oPj = 0- Therefore, even-parity methods are not zero-stable and we can 
restrict our attention to odd-parity methods, c = — 1. 



If the reversibility criterion (|35| ) is satisfied, then 

L[x{t), h] = cL[x{t + kh), -h]. (37) 
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n=0 "'■ 




m 


We can re-write this equation as 






n=l ''" 




(39) 


If c(— 1)*' = —1, then this equation imphes that Cp = ii Cq = ■ ■ 


■ = Cp., 


= and thus the 



order cannot be p — 1. Therefore, odd-parity methods have even order. 

The stabihty of multistep methods for oscillatory problems can be parametrized by the 
interval of periodicity introduced by Lambert & Watson (1976). Suppose the right-hand size 
of equation (|I]) is f{x,t) = iujx, with uj real. Then, the multistep method (|2l|) becomes a 
linear difference equation, with solution Xn = a^"- Here, a and i^ are complex constants, the 
latter satisfying 

p{0 - ^ooha{0 = or g(e) ^ -t^ = ooh. (40) 

(7[e ) 

The interval of periodicity is the largest value of uh such that all of the roots of the first of 
equations ( |iO| ) lie on the unit circle. Outside the interval of periodicity, the solution grows 
exponentially and hence is unstable. For reversible methods with odd parity, we may write 

... _ E,to«,sin[(j-lfc)g] 
'^^ E,to/5.cos[(j-ifc)^]- ^''^ 

If the multistep method is zero-stable, then p(^) has k distinct roots on the unit circle. 
Thus, g{d) has k distinct roots on the interval [— 7r,7r]. For sufficiently small values of uh, 
the equation g{6) = ujh will still have k roots. As uh is increased, eventually a pair of these 
roots will disappear. This occurs at the smallest local maximum of g{9) and marks the end 
of the interval of periodicity. So, a plot of g{6) can be used to determine the interval of 
periodicity (Fukushima 1998). 



2.2. Instabilities in reversible multisteps 

The growth of errors in multistep methods is discussed by Henrici (1962), Gear (1971) 
and especially by Cano & Sanz-Serna (1998). For our purposes, a limited heuristic version 
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of these analyses is sufficient. Let {xj} be a solution of the multistep method (pi]). Now 
perturb the solution to {xj + Cj}, \ej\ ^ \xj\. Linearizing equation (pT]), we have 



J2{ak-j - hpk-jfL-j)en~j = 0, 



(42) 



i=o 



where /' = df{xj)/dx. Now look for solutions of the form e„ = ^^y{to + nh) where y{t) 
is smooth and ^ is a complex constant. In the limit h ^ 0, Xn —>■ x{tQ + nh), where x{t) 
is the accurate solution of the differential equation (|l]) and therefore is smooth. Then ( ^21) 
becomes 

' dx 



k 

E 

i=o 



aj — h[3j^{x{t + jh)) 



ey{t + jh) = o, 



(43) 



where t = tQ + {n — k)h. Now expand in a Taylor series in h. To order h^ we have y{t)p{C,) = 0, 
which requires that ^ is one of the roots ^j,j = 1, . . . ,k oi the polynomial p. To order h we 
have 



dyjt) 
dt 



J2j=oPj<.^ \ u\'^f l~U\\ 



x,=oJ«.e^ 



9x 



^(0 
ep'(0 



|(^W). 



(44) 



Thus, as /i — i> 0, the linearized difference equation ( ^2]) has solutions of the form 






(45) 



where yj{t) satisfies the differential equation 

dyj _ , df{x{t)) 



dt 



Ar 



dx 



Vjii), 



and the growth parameters are 



A, 






3 



k. 



(46) 



(47) 



This should be contrasted with the variational equation of the orbit itself. If the trajectory 
x{t) is perturbed to x{t) + e(t), \e\ <^ \x\, then 



de{t) df{x{t)\~^^~^ 



dt 



dx 



(4J 



Cano & Sanz-Serna (1998) point out that the variational equations (^6]) and ( ^8]) generally 
have quite different solutions, and hence the former can be unstable even if the latter is 
stable. However, this instability can be evaded for special values of the growth parameters 
1, equations (0) and ( ^81) are the same, for Aj = —1, equation ( ^Tj) is simply 



Xf for Aj 
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the time-reversed version of equation (^81), and for A^ = 0, the solution of equation (0) is 
simply yj =const. We use the term "zero-growth" methods to denote multistep methods 
such that \j G { — 1, 0, 1}, j = 1, . . . , /c. Cano (1996) has shown that for the reversible planar 
Kepler problem, these are the only possible choices for the growth parameters \j to ensure 
stability, although for other potentials this may not be the case. 

Multistep methods for integrating the special second-order equation x" = f{x,t) (Lam- 
bert & Watson 1976; Quinlan & Tremaine 1990) are not generally subject to this sort of 
instability, for the following reason. The expansion of the analog to equation (^31) in this case 
yields yit)p{^) = to order h^ and y'(t)p'{^) = to order h. Satisfying these two equations 
simultaneously requires that ^ is a double root of p{C)', therefore there is no instability if all 
of the roots of p(^) are simple (Cano & Sanz-Serna 1998). 



3. Zero-growth methods 

In this section, all the zero-growth multistep algorithms with up to six steps are found. 
Let us begin with some general results for stable, odd-parity, multistep RiAs, specified by 
the characteristic polynomials p(^) and cr(^). Let ^j,j = 1, . . . , A; be the roots of p(Oj ^oi 
zero-stable reversible methods, these are distinct, lie on the unit circle and are either real 
(+1 or —1) or appear in complex-conjugate pairs. Equation ( ^9]) implies that 1 is always 
a root, which we denote ^i. Therefore, —1 is also a root if and only if the total number of 
roots k is even; when present, we denote this root as .^2- We normalize the multistep method 
so that Ok = 1 and 

k 

Pio = m^-Q- (49) 

We now demand that the growth parameter associated with the root ^i is A^ G { — 1,0,1} 
(note that Ai = +1 by eq. |5D|). From equation (^Tj), this requires that 

^{^i) = >^i^ip'{^i) = >^Al[(^i-Q- (50) 

The unique polynomial of order k that passes through the points (^) and satisfies the 
symmetry condition (^) is 

k 

^(0 = iEAKe + eon(e-0), (51) 

1=1 j^i 

with Xj = Xi when ^j^i = 1. An alternative form is 

^(0 = MOEA.|^- (52) 
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Note that if A^ = then ^^ is a root of both p(^) and cr(^). In this case, there is a 
zero-growth multistep RIA with the same order and error constant but fewer steps (see §0). 
Therefore, we restrict our attention to methods with A; G {—1, !}• 

The method is exphcit if a{0) = 0, which in turn requires 

k 

E ^J = 0' (53) 

a condition that can only be satisfied if k is even. 

A /c-step method is completely specified by the fc*'^-order polynomials p{C,), <^{0- ^^ ^ 
stable multistep RiA with k even, the roots of p(^) are +1, —1, Xj ± ii/j where x'j + y'j = 1, 
j = 1, . . . |(/c — 2). For k odd, the root at —1 is not present and the index runs to |(A; — 1). 
In a zero-growth method, the polynomial o"(^) is determined by equation (^) once p(^) and 
the signs Xj are specified. We now discuss all stable, zero-growth, multistep RiAs of < 6 
steps. We label the interesting methods by SZ/cp where k is the number of steps and "p" is 
an optional suffix that distinguishes different methods with the same k. 



3.1. One- and two-step methods 

Since ^i = 1 is always a root, for /c = 1 we must have p(^) = ^ — 1. Since Ai = +1 by 
equation (|30|) , the zero-growth requirement is automatically satisfied. Equation (^) requires 
that cr((^) = 1(^ + 1), which yields the trapezoidal method (eq. |T5|) 

Xn+i =Xn + \h{fn+i + /„); method SZl (54) 

The first few terms of the Taylor series ( |2^ are 

Co = 0, Ci=0, C2 = 0, C3 = -j^; (55) 

thus the error constant is C = — ^ and the method is second-order. The function g{9) (eq. 
40D is 2 tan |6', so the interval of periodicity is infinite. 



For k = 2 the roots of p(^) can only be ±1 (^i = 1 is always a root; ^2 niust be real 
since ^2 is also a root and there are only two roots; ^2 must be on the unit circle and distinct 
from ^1; thus ^2 = ~1)- Thus, the first characteristic polynomial is 

p(o = e-i. (56) 

There are two choices for the second characteristic polynomial cr(0) depending on A2 G 

{-1,1}: 

a(0 = 2e, or a{0 = e + l- (57) 
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Fig. 1. — The fractional energy error after integrating the standard Kepler problem (a = 1, 
GM = 1) for 100 time units using the 2-step RIA (|6T|) , as a function of the parameter Pq. The 
two solid curves were computed with timesteps h = 0.01 and 0.001 for eccentricity e = 0.2, 
while the dashed curve has h = 0.001 and e = 0.5. Off-scale energy errors indicate that an 
implicit timestep did not converge after 20 iterations. The minima at /3o = 0, |, 1 correspond 
to the explicit midpoint method, the trapezoidal method, and the trapezoidal method with 
a timestep of 2h. 
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The first of tliese gives tlie explicit midpoint metliod (eq. 0), 

Xn+i = Xn-1 + 2/i/„ metliod SZ2. (5^ 



The first few terms of the Taylor series (|2^ are 



Co = 0, Ci = 0, C2 = 0, ^3 = 1; (59) 

thus the error constant is C = | and the method is second-order. The function g{6) is sin 6' 
so the interval of periodicity is (ti;/i)max = 1- The second equation for cr(^) in ( ^7|) yields 



Xn+l=Xn-l + h{fn+l + fn-l), (60) 

which is the same as the trapezoidal method (|5^ , with a timestep of 2h. 

All of these two-step methods are special cases of a one-parameter family of multistep 
RIAS, 

Xn+l = Xn-\ + /i[/5o/n+l + 2(1 - /3o)/n + /^o/n-l]- (61) 

The explicit midpoint method (eq. |SS|) has /^o = and the trapezoidal method with timestep 
2h (eq. |6^) has /5o = 1. When /?o = |; the method is simply the composition of two 
trapezoidal steps (eq. ^^. Milne's method (/5o = |) is the two-step method of maximum 
possible order (p = 4), but its disadvantage is that its growth parameter A2 = — | and so it 
is susceptible to numerical instability. The interval of periodicity (co'/i)max = (1 — 2/3o)^^^^ 
for /5o < I and infinity for /5o > |. To compare these methods we have integrated a Kepler 
orbit with unit semimajor axis and period 27r for 100 time units using various eccentricities 
and timesteps (Fig. |1|). Even for timesteps as low as /i = 0.001, the relative energy error 
\IS.E\IE ^ 1, except near the zero-growth methods at Pq = 0,^,1. In particular, the 
second-order zero-growth methods - even the explicit method at /?o = - exhibit much 
better behavior than the fourth-order Milne's method. 



3.2. Three- and four-step methods 

For k = 3, the distinct roots of ^ can only be ^1 = 1 and ,^2,3 = u±iv where \u\ < 1 and 
v = (1-m2)V2. Thus, 

p(0 = e- (2w + l)e + i2u + l)e - 1. (62) 

There are two choices for o"(^). The first choice is to take A2 = A3 = -|-1, so that we have 

a(0 = le - (I + u)e - (i + «)^ + §• (63) 
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The first few terms of the Taylor series (E^ are 



Co = 0, Ci = 0, ^2 = 0, Cs = -f + lu. (64) 

C3 = has no solution for \u\ < 1 so the method is second-order; the interval of periodicity 
is infinite. The second choice is to take A2 = A3 = — 1 so that we have 

a(o = ^le + (I - u)e + (I - u)^ - \. (65) 

The first few terms of the Taylor series (^4]) are 

Co = 0, Ci = 0, C2 = 0, C3 = ^ + |m. (66) 

Once again, C3 = has no solution for \u\ < 1 so the method is second-order. The interval 
of periodicity ranges from 2~^/^ as m — i> — 1 to as m — ;> 1. 

For A; = 4, the distinct roots of ^ can only be ^1 = 1, ^2 = —1 and ^3,4 = u ±iv where 
\u\ < 1 and f = (1 - M^)^/^. Thus, 

p(0=e'-2<3 + 2<-l. (67) 

There are four choices for o"(^): A2 G { — 1,+1}, A3 = A4 G { — 1,+1}. None of these yield 
methods of order > 2 for |u| < 1. 

The three- and four-step methods have no obvious advantage over the trapezoidal or 
explicit midpoint methods, and we will not explore them further. 



3.3. Five- and six-step methods 

For k = 5, the distinct roots of ^ are ^1 = 1, ,^2,3 = mi ± ivi and ^4^5 = M2 ± iv2, where 
\uj\ < 1 and Vj = (l — u'^Y^'^. There are three choices for a{^): (A2, A4) = (—1, —1), (—1, +1), 
or (+1, +1) (note that (+1,-1) is not a distinct choice), with A3 = A2 and A5 = A4. The 
choices (A2, A4) = (—1, —1), (+1, +1) yield only second-order methods. The choice (A2, A4) = 
(—1, +1) yields a fourth-order method if 



U2= ,, ^. , «iG(-l,l), M2G(-f,l). (68) 

The integration formula is 



Xn+i = (x„ -a;„_3)(l + 2ni + 2^2) -2(a;„_i -x„_2)(l + ^ii + ^2 + 2^1^2) + a;„_4 
+ |/i[/„+i + (/„ + /„_3)(1 + 2ui - 6U2) 
+2(/n-i + /n-2)(l-3Mi+M2 + 2MiM2) + /n-4] method SZ5. (69) 
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The error constant is 

(t(1) 1440(mi-1) ^ ^ 

For k = 6, the distinct roots of ^ are ^i = 1, ^2 = — 1, ^3,4 = ui Hvi and ,^5,6 = U2±iv2, 
where \uj\ < 1 and Vj = (1 — u'^Y^'^. There are six choices for cr(^): A2 = —1 or +1 and 
(A3, A5) = (-1, -1), (+1, -1), (+1,+1), with A4 = A3 and Ae = A5. The choices (A3, A5) = 
(—1, —1) or (+1, +1) yield only second-order methods. The choice (A2, A3, A5) = (+1, +1, —1) 
yields a fourth-order method if 

M2 = ^^^, Ml G (-1,1), M2G(-i,l). (71) 

4 — Ml ^ 

The error constant is 

C=^= ^f + "\. (72) 

cr(l) 45(mi-1) ^ ' 

The integration formula is 

Xn+i = 2(xn- Xn-4) (mi + U2) - (a;„_i - X„_3) (1 + 4M1M2) + a^n-5 

+ h[fn+l + fn-5 - 4(/„ + fn-i)u2 

+ (/„_i + /„-3)(3 + AU1U2) - 8fn-2Ui] method SZ6i, (73) 

the label "i" stands for "implicit" . 

The choice (A2, A3, A5) = (—1, +1, —1) yields an explicit fourth-order method if 

7mi-1 
ui + 5 

The error constant is 



"2 = ^:^^^, Kie(-i,l), K2e(-i.i). (74) 



q ^ 19 + 11.. 

a{l) 180(1 -Ml) ^ ^ 

The integration formula is 

Xn+i = 2(x„-x„_4)(mi + M2) - (a;„_i -x„_3)(l +4M1U2) +a;„_5 

+ h[2{fn + /„-4)(l + Ml - M2) - 4(/,_i + /„_3)(Mi + M2) + 

4/„_2(l — Ml + M2 + 2M1M2)] method SZ6e. (76) 

The intervals of periodicity for the three fourth-order methods are shown in Figure ^ These 
are largest for mi near —1 and decrease rapidly for mi > 0. Thus, only methods with negative 
Ml are of practical use. At best, the interval of periodicity is significantly smaller than that 
of comparable multistep RiAs for special second-order equations (Lambert & Watson 1976; 
Quinlan & Tremaine 1990; Fukushima 1998). This is the price of the greater generality of 
the present methods. 
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Fig. 2. — The interval of periodicity for three zero-growth, fourth-order RiAs: the five-step 
method SZ5 defined by equations ( |68D and (|69D , the 6-step imphcit method SZ6i defined by 
(fn]) and (|73D , and the 6-step exphcit method SZ6e defined by (0) and (|76D . 
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Fig. 3. — The largest timestep that can integrate a circular orbit in the reversible Kepler 
potential and maintain fractional energy error below 10~^ for 100 orbital periods. Results 
are shown for the three zero-growth, fourth-order RiAs: the five-step method SZ5 defined 
by equations ( |BBD and (p^), the 6-step implicit method SZ6i defined by ( [7TD and ([75|), and 
the 6-step explicit method SZ6e defined by (|7^ and (|76|). Although results are displayed for 
circular orbits, the curves for mildly eccentric orbits (e ^ 0.1) are almost indistinguishable. 
This diagram can be compared with Figure ^, which shows the theoretical upper limit to 
the timestep for the harmonic oscillator potential. 
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4. Numerical experiments 

This section supplements our theoretical analysis of multistep RiAs with numerical ex- 
periments using the Keplerian Hamiltonian |f ^ — 1/r. The equations of motion are 

where r^ = x^ + y^ and the gravitational constant G and the attracting mass M have been set 
to unity. Our test integrations follow an orbit with eccentricity e, unit semi-major axis, and 
orbital period 2tt. The orbit is started at apocenter on the x-axis, so the initial conditions 
are 

x = l + e, i/ = 0, t;, = 0, Vy = [{l-e)/{l + e)Y'\ (78) 

The interval of periodicity shown in Figure ^ gives the maximum stable timestep for the 
harmonic oscillator potential only. It is of greater relevance to astronomers to establish the 
timesteps that can be safely used to integrate circular or mildly eccentric orbits in the Kepler 
potential. Figure ^ shows the largest timestep for which the relative energy error is less than 
10^^ after 100 orbital periods. The results displayed are for circular orbits, but the curves for 
mildly eccentric orbits (e < 0.1) are almost indistinguishable. A comparison between Figures 
and reveals that the timestep needed for accurate integration of Keplerian near-circular 
orbits is typically at least a factor of five smaller than would be naively inferred from the 
interval of periodicity. 

Figure ^ shows a further test of the multistep RiAs. The standard Kepler problem is 
integrated for 100 time units with timestep h = 0.003 and eccentricities 0.2 and 0.5. A 
striking feature is that the error is almost independent of the parameter ui over a limited 
range, while outside this range the method is unstable. The stable range shrinks as the 
eccentricity increases. The dependence of energy error on timestep is explored in Figure ^ 
In this Figure, we have specialized to a single value of the parameter ui. —0.75 for SZ5 and 
SZ6i and —0.25 for SZ6e. All three methods display the expected dependence AE oc h'^ 
for h < 0.01. For larger values of h, the methods are unstable. We have also plotted 
the behavior of the classic fourth-order Adams-Bashforth (explicit) and Adams-Moulton 
(implicit) multistep methods. The classic methods exhibit similar errors for h < 0.01 but 
remain stable at much larger timesteps. The advantage of the reversible multistep methods 
only becomes clear over longer time intervals. Figure ^ shows that the maximum energy error 
in the reversible methods is constant at large times, while the energy error in the classical 
methods grows linearly. 
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Fig. 4. — The fractional energy error after integrating the standard Kepler problem (a = 1, 
GM = 1) for 100 time units using timestep h = 0.003 and eccentricities 0.2 (left panel) 
and 0.5 (right panel). The figure shows results for three fourth-order methods: the five-step 
method SZ5 (short-dashed lines); the 6-step implicit method SZ6i (long-dashed lines), and 
the 6-step explicit method SZ6e (solid lines). The curves for e = 0.5 are always higher than 
the corresponding curves for e = 0.2. 
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Fig. 5. — The fractional energy error after integrating the standard Kepler problem (a = 1, 
GM = 1) with eccentricity e = 0.2 for 100 time units. The figure shows results for three zero- 
growth, fourth-order RiAs: the five-step method SZ5 with ui = —0.75 (short-dashed line); 
the 6-step implicit method SZ6i with ui = —0.75 (long-dashed line); and the 6-step explicit 
method SZ6e with ui = —0.25 (solid line). Errors are also shown for two other fourth-order 
multistep methods: Adams-Bashforth (dot-short dashed line) and Adams-Moulton (dot-long 
dashed line). 
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Fig. 6. — The maximum fractional energy error over integration time t for the standard 
Kepler problem (a = 1, GM = 1) with eccentricity e = 0.2 and timestep h = 0.005. The 
figure shows results for three zero-growth, fourth-order RIAs: the five-step method SZ5 with 
Ml = —0.75 (short-dashed hue); the 6-step imphcit method SZ6i with ui = —0.75 (long- 
dashed line); and the 6-step explicit method SZ6e with ui = —0.25 (solid line). Errors 
are also shown for the fourth-order Adams-Bashforth (dot-short dashed line) and Adams- 
Moulton methods (dot-long dashed line). 
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4.1. Variable timesteps 

Variable timesteps are useful for eccentric Kepler orbits. To incorporate these, we 
introduce a fictitious time r by the relation dt = g{x,y)dT, so the equations of motion 
become 

dx dy dvx x dvy / \ V dt 

^ = ^(x,i/K j^ = gix,y)vy, -jp = -9i^^y)-s^ -d7 = -'^'^y^7^^ ^7 = ^^^'^^- 

(79) 

These are integrated using unit timestep in r. We typically use g{x,y) oc r^/^ since this is 
the characteristic free-fall time from radius r. We parameterize the timestep by the number 
of force evaluations per unit time (or per radian, since the orbital period is 2tt). This 
is distinct from the number of steps per radian because implicit methods require several 
iterations per step, and because some explicit methods - though not multistep methods 
- require several force evaluations per step. Figure |^ shows the fractional energy errors 
resulting from integrating Kepler orbits using the trapezoidal method (eq. ^ solid lines) 
and explicit adaptive leapfrog (eq. |1^, dashed lines). Each step of the trapezoidal method 
was iterated to convergence, with the first approximation taken from the first-order Euler 
method - typically 3 iterations were required at the smallest timesteps, rising to ~ 10 at 
the largest. The energy errors are the maximum over 1000 time units or SOO/vr orbits. The 
results for shorter integrations over 100 time units are almost indistinguishable, showing that 
there is no secular energy drift with either integration method. Adaptive leapfrog generally 
gives energy errors that are smaller by about an order of magnitude, mostly because it 
is explicit so there are fewer force evaluations per timestep. We have conducted similar 
experiments with the explicit midpoint method (eq. ^ not shown in Figure) but this is 
much less successful at following high-eccentricity orbits, presumably because of its smaller 
interval of periodicity. 

We next investigate the fourth-order multistep RiAs. These have relatively small in- 
tervals of periodicity and cannot reliably integrate high-eccentricity orbits. However, they 
are successful at following orbits with moderate eccentricities. Figure ^ shows the fractional 
energy errors resulting from integrating a Kepler orbit with eccentricity e = 0.5 using four 
methods: the explicit six-step method with ui = —0.25 (method SZ6e, eq. |7^); the implicit 
six-step method with ui = —0.75 (method SZ6i, eq. |73D; the implicit five-step method with 
Ml = —0.75 (method SZ5, eq. ^^. In all panels, the heavy solid line is the error resulting 
from integrating the eccentric orbit with variable timestep, the light solid line is the error 
from integrating the same orbit with fixed timestep, and the dashed line is the error from 
integrating a circular orbit with fixed timestep. For large timesteps, the implicit methods 
do not converge and in this case no energy error is plotted. For all four methods a variable 
timestep yields errors that are 1-2 orders of magnitude smaller than a fixed timestep (at the 
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Fig. 7. — Results of integrating the differential equations for a Kepler orbit (ffOD for 1000 
time units, using the trapezoidal method (eq. ^, solid lines) and explicit adaptive leapfrog 
(eq. |1^, dashed lines). The orbit eccentricities are e = 0.5 (stars), 0.9 (triangles), 0.99 
(squares), 0.999 (pentagons), 0.9999 (circles). The vertical axis is the maximum fractional 
energy error during the integration and the horizontal axis is the number of evaluations of 
the right-hand side of equations (|79|) per unit time (the orbital period is 27r). 
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Fig. 8. — Integration of a Kepler orbit with eccentricity e = 0.5 by four fourth-order methods: 
the exphcit six-step method with ui = —0.25 (method SZ6e, eq. ^); the imphcit six- 
step method with ui = —0.75 (method SZ6i, eq. |7^); the imphcit five-step method with 
Ml = —0.75 (method SZ5, eq. B9|); and the imphcit, fourth-order, two-stage Runge-Kutta 



RIA defined by equations (TT). The heavy sohd fines are the fractional energy errors for 
e = 0.5 and variable timestep, the light solid lines are for e = 0.5 and fixed timestep, and 
the dashed lines are for e = and fixed timestep. 
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same number of function evaluations per unit time). The explicit multistep method performs 
much better than the implicit methods, in part because the implicit methods require 3-8 
force evaluations per step to iterate to convergence. However, even if the convergence of 
the implicit multistep methods could be achieved in 2-3 iterations they would not perform 
significantly better than the explicit method on this test. A surprising result is that there is a 
range of the horizontal axis in which the multistep methods with variable timestep integrate 
eccentric orbits with smaller energy error than circular orbits. 



4.2. Resonant instabilities 

Toomre (1990, private communication) and Quinlan (1999) showed that multistep RiAs 
for the special second-order differential equation (|19D suffer instabilities at unlucky timesteps 
at which there is a resonance between the oscillation frequencies of the solution and the 
method. This phenomenon must also occur in the multistep RiAs for first-order differential 
equations that are investigated here. Suppose the angular speed of the circular orbit given 
by the numerical method is tu. If the number of steps in the multistep RIA is at least 5 
(or 6 for second-order differential equations), then there are two or more complex-conjugate 
pairs of roots of the characteristic polynomial (|2^) on the unit circle. These we write as 
^i = exp{±i6i) and C,j = exp{±i6j). Quinlan (1999) shows that the troublesome timesteps h 
satisfy 

Oi - 9j = 2uh. (80) 

In practice these resonant instabilities are never a concern for our zero-growth multistep 
RiAs. This is because the timesteps (|80D lie outside the interval of periodicity for nearly 
circular Keplerian orbits. Let us demonstrate this by considering one example in more 
detail - namely, the 6-step implicit method SZ6i defined by ( fflD and (^). This is the 
most susceptible to resonant instabilities of the zero-growth methods we have examined, 
both because it has the largest interval of periodicity and because, as Ui —>■ —1, two of 
the spurious roots coalesce as ui -^ —1. We recall that there are five spurious roots on 
the unit circle, namely ^2 = —1, and the two complex-conjugate pairs ^3^4 = mi ± ivi and 
^5,6 = M2 ± 'i'V2, where Vj = (1 — u'^Y^'^. Using ([7I|), the angular separations between the 
spurious roots on the unit circle can be readily computed as a function of Ui. Figure ^ 
shows the locations of the potentially troublesome timesteps h, together with the maximum 
possible timestep (already presented in Figure ^. The resonant instabilities always occur 
at timesteps that are outside the interval of periodicity for Keplerian orbits. Similar results 
hold for the methods SZ5 and SZ6e. 
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Fig. 9. — For the 6-step implicit method SZ6i, the resonant timesteps are shown in the full 
lines. The interval of periodicity for nearly circular Keplerian orbits is shown in a dashed line. 
The resonant timesteps are always larger than the maximum permitted timestep. Similar 
results hold true for SZ5 and SZ6e. 
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5. Conclusions 

The generic instabilities of multistep RiAs for first-order differential equations identi- 
fied by Cano & Sanz-Serna (1998) are a grave problem and render many of these methods 
unusable. The main contribution of this paper is the elaboration of a class of multistep 
RiAs - the zero-growth methods - that successfully evade the instabilities. This select group 
includes two familiar second-order multisteps, namely the trapezoidal method (|5^ ) and the 
explicit midpoint method (0), as well as entirely new multistep RiAs. In particular, we have 
identified three one-parameter families of zero-growth, fourth-order RiAs - namely, the five- 
step implicit method SZ5 defined by equations ( |68D and (p9|), the six-step implicit method 
SZ6i defined by (|7TD and (|73|), and the six-step explicit method SZ6e defined by ([7^) and 



([761). Of these, the explicit method SZ6e eliminates the need to iterate to convergence at 
each timestep without a major sacrifice in the interval of periodicity. The stability of the 
zero-growth methods has been understood at a theoretical level and confirmed at a practical 
level by integrations of circular and eccentric Keplerian orbits. 

There are several advantages of multistep RiAs for first-order differential equations over 
the multistep RiAs for special second-order differential equations examined by Lambert & 
Watson (1976), Quinlan & Tremaine (1990), and Fukushima (1998, 1999). They permit the 
introduction of variable timesteps without spoiling reversibility or long-term energy conser- 
vation. They are also more general, since any ordinary differential equations can always be 
reduced to a system of first-order equations; thus, for example, the methods described here 
can be used to follow orbits in rotating frames of reference or the motions of rotating rigid 
bodies. A third advantage is that they are less susceptible to timestep resonances. The price 
paid for these attractive features is that (i) the interval of periodicity of the higher order 
methods such as SZ5, SZ6i and SZ6e is rather small, forcing the use of smaller timesteps than 
other methods; (ii) more steps are needed for a method of given order; there are no fourth- 
order methods with fewer than five steps and no fifth or sixth-order methods with 8 or fewer 
steps (the limit of our explorations). Nonetheless, the six-step explicit method (SZ6e) re- 
mains a competitive option for problems in which it is critical to have variable timesteps and 
to avoid irreversible numerical errors: examples include long-term integrations of asteroid or 
comet orbits, or orbits near the centers of galaxies. 

There are several possible avenues for future exploration: 

1. Are there zero-growth RiAs with more than k = 6 steps that have higher order or a 
larger interval of periodicity? We have explored all k < 8 but have found no RiAs of 
order greater than 4. 

2. Gautschi (1961) developed a class of methods closely related to the linear multisteps. 
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Just as a linear multistep is defined by the requirement that the operator ( ]23|) anni- 
hilate all algebraic polynomials of order < p, so the Gautschi multisteps annihilate 
certain trigonometric polynomials up to a given degree. These methods are particu- 
larly appealing when the solution is known to be periodic and a reasonable estimate 
for the period of the orbit can be guessed in advance - requirements that are often 
satisfied in solar system integrations. Can we find reversible Gautschi multisteps? 

Integration methods are generally designed to maximize the order p, defined so that 
the one-step error in following a system with characteristic frequency u: is proportional 
to {hujy~^^. Perhaps it is more sensible to "economize" the error, by minimizing the 
maximum value of the error over a range of frequencies < u; < cUmax- In this case the 
error would be nearly proportional to a Chebyshev polynomial Tp+i(x) with argument 
X = tu/coimax and degree p + 1, instead of a;^"*"^. 
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